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FIRST-AND SECOND-ORDER SENSITIVITY ANALYSIS OF A P- VERSION FINITE ELEMENT 
EQUATION VIA AUTOMATIC DIFFERENTIATION 


Abstract 

Sensitivity analysis is a technique for determining derivatives of system responses with respect to design param- 
eters. Among many methods available for sensitivity analysis, automatic differentiation has been proven through 
many applications in fluid dynamics and structural mechanics to be an accurate and easy method for obtaining deriv- 
atives. Nevertheless, the method can be computational expensive and can require a high memory space. This project 
will apply an automatic differentiation tool, ADIFOR, to a p - version finite element code to obtain first- and second- 
order thermal derivatives, respectively. The focus of the study is on the implementation process and the performance 
of the ADIFOR-enhanced codes for sensitivity analysis in terms of memory requirement, computational efficiency, 
and accuracy. 


1. Introduction 


Response derivatives are important to many engineering applications, such as design approximation, design pre- 
diction, design optimization, etc. Methods have been developed in the past to calculate response derivatives systemat- 
ically. Finite differencing, direct differentiation, semi-analytical method, and adjoint variable approach are a few 
examples that are applicable to the algebraic system equations. Users of those methods are required to generate com- 
puter codes according to the sensitivity equations derived by using calculus. Automatic Differentiation (AD), on the 
other hand, works directly with the computer compiler to generate an enhanced code that includes the derivative 
information of the responses computed in the input computer code 1,2 . 

Quite a few publications have documented the experience of using AD tools for sensitivity applications. One AD 
tool, ADIFOR (Automatic Differentiation of FORtran) 3 , has been used to differentiate a finite element structural code 
to obtain first-order derivatives 4 . ADIFOR has also been applied to commercially rated computational fluid dynamics 
codes to obtain first-order aerodynamic sensitivities 5,6 . Particularly in Ref.7, ADIFOR is not directly applied to the 
entire computational Fluid Dynamics (CFD) code as a black box; instead, it is applied to the selective subroutines that 
define the residuals of the aerodynamic algebraic equations. The same procedure is extended to find the second-order 
aerodynamic derivatives 7 . Another application of ADIFOR can be found in Ref.8 for flutter-speed sensitivity analy- 
sis. Valuable observations are made in the paper regarding the use of ADIFOR for sensitivity analysis. These 
publications 4 ' 8 have illustrated that the ADIFOR-enhanced codes are accurate for sensitivity analysis and can be gen- 
erated with minimal coding effort. Nevertheless, they observed that the ADIFOR-enhanced codes can be bulky and, 
compared to the hand-coded sensitivity module, may require more computer resources to compute and store the sen- 
sitivity information. 

The current work will apply ADIFOR to a ^-version finite element code for one-dimensional heat transfer. The 
study will investigate the implementation process, the memory requirement of the ADIFOR-enhanced codes and the 
computational efficiency and accuracy of the overall sensitivity analysis using the ADIFOR-enhanced codes. The 
focus is placed particularly upon the second-order sensitivity analysis, which is known for its complexity in deriva- 
tion and its intensity in computation. 


2. Problem Formulation 


The code of concern is a p - version finite element which was written by Bey and Wilson 9 . The code is designed to 
solve one-dimensional, steady-state heat transfer problems. Besides the conventional constant and linear polynomi- 
als, high order polynomials, with zero values and derivatives at the end nodes, are chosen as interpolating functions . 
In this way, the domain to be analyzed can be discretized in the conventional way, where the nodal temperatures are 
associated with the constant and linear polynomials. The accuracy of the solution can be improved by increasing the 
number of “nodeless” degrees of freedom associated with the higher order polynomials. 

The finite element equations for the thermal codes can be expressed symbolically as 
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Where K, «, and/ in represent the heat conduction matrix, the temperature vector, and the vector of external heat gen- 
eration, respectively. The dimension of the heat conduction matrix K depends upon the order of the polynomials used 
in the model. 

The first-order sensitivity equation can be obtained by differentiating Eq. (1) with respect to an arbitrary design 
parameter, b it to obtain 
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The process can be continued to obtain the second-order sensitivity equation 
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where bj is another design parameter. 

Equations (2) and (3) provide the values of au/36, and d 2 u/db i db j for computing derivatives of an arbitrary 
performance function, t|/(w) , required in the gradient-based design optimization process: 
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where the superscript T represents the transpose. Those derivatives of V(m) can also be calculated by the adjoint 
variable method 14 as 


and 
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where the first-order derivative, 3 u/dbj is obtained by Eq. (2) and the adjoint variable vector, A. , is obtained by solv- 
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ing the adjoint variable equation 



( 8 ) 


Note that the sensitivity equations, Eqs. (2) and (3), and the adjoint variable equation, Eq. (8) can be solved with 
backward substitution, as matrix K has been factorized already when Eq. (1) is solved. Nevertheless, solving sensi- 
tivity equations remains a difficult computational task. Many new derivative terms not found in analysis are needed 
now to construct either the right-hand sides of the sensitivity equations, Eqs. (2) and (3), or the derivatives of Eqs. (6) 
and (7). These new terms include the derivative dK/db a three-dimensional matrix, and the derivative 5 K/dbfib ■ , 
a four-dimensional matrix. Furthermore, the number of sensitivity equations to be solved can be large. In fact, m and 
( m{m + 1 )/2 ) numbers of linear equations are needed respectively to solve for the first- and second-order deriva- 
tives of u or \}/(«) in Eqs. (2)-(5), where m is the number of the design parameters. As for the adjoint variable 
method, n and (m + n) numbers of linear equations are needed in Eqs. (6)-(8), where n is the number of the perfor- 
mance functions. Although using the adjoint variable method may reduce the number of linear equations solved for 
sensitivity analysis in the case where n is less than m , computing the new derivative terms in the sensitivity equa- 
tions is still required. 

Though the design parameters, b t , are directly related to the stiffness matrix, and may also be dependent on each 
other. In fact, they may be “linked” to the true design variables, x , through a pre-determined relation, b(x ) . In this 
case, for a specific design variable, x i , Eqs. (2) and (3) are rewritten as 
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where s, denotes db/ dx., which is often called the seed matrix in the ADIFOR literature and is usually provided by 
the user. Equations (6) and (7) can be modified in a similar manner to account for the design variable linkage. 

3. Implementation Procedure 


The main task of implementation is to use ADIFOR to construct the new derivative terms in the sensitivity equa- 
tions. ADIFOR is used here selectively to minimize the coding effort and to minimize the memory requirement while 
maximizing the computational efficiency of the ADIFOR-enhanced code. The step-by-step procedure is described 
below for the first-and second-order sensitivity analyses. The procedure is the same for every design variable; hence, 
explanation of each step is given only for the shape design variables. 

The first step involves the collection of the FORTRAN subroutines that compute the elemental stiffness matrix, 
K e . These subroutines are input to ADIFOR to obtain the derivatives of the element stiffness matrix with respect to 
the nodal coordinates, b e , of a typical element, which are specified as independent variables. The output of ADI- 
FOR gives a collection of new subroutines that can compute K e as well as (( dK e / db e ) S ) where S is a user-speci- 
fied seed matrix. In this step, S is fixed as an identity matrix with which the ADIFOR-produced code computes 
dK e /db', as a three- dimensional array. 

The second step generates a new subroutine to compute the product (dK e /db e )u e where u c is the solution of 
Eq. (1) pertaining to the element under consideration. The product is a two-dimensional array. 
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The third step generates another new subroutine to compute the product (( dK e /db e )u e )S e where S c is a user 
specified seed matrix, db/dx , derivatives of nodal coordinates with respect to independent geometric design vari- 
ables. 

The fourth step assembles the element arrays produced in the third step into global arrays which represent the 
first term on the right-hand side of Eq. (9). 



Fig.l Computational procedure for analysis and 
sensitivity analysis. 

The second term on the right-hand side of Eq.(9) and the second and the third terms of Eq. (10) can be con- 
structed in a similar manner. The above procedure, however, needs further work in order to compute the most compli- 
cated term: the first term on the right-hand side of Eq. (10). One may proceed to submit the code^produced in step one 
to ADIFOR once more to obtain the second- order derivatives with a seed matrix, ( d K e / db e )S e . Again the seed 
matrix is specified as an identity matrix. The u e is then multiplied to the second-order derivatives as done in the above 
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Step two. The resultant three dimensional arrays are then multiplied by S e twice to obtain 
((( d 2 K e /db 2 e )u e )S e )S e where S e is db/dx . The above fourth step is then followed to complete the first term on 
the right-hand side of Eq. (10). The code to compute the fourth term in Eq. (10) can be generated in a similar manner. 

Once the computer codes that compute the right-hand sides of Eqs. (9) and (10) are assembled, they can then be 
attached to the original analysis code as a separate module for sensitivity analysis. Such an arrangement is depicted in 
the flow chart of Fig. 1 . 

To study the performance of ADIFOR, numerical experiment is done to compare memory requirement and com- 
putational efficiency of using the ADIFOR-produced code versus the finite difference method for sensitivity analysis. 
The results are listed in Table 1 for the first- and second-order shape derivatives of a one-dimensional, p- version heat 
conduction element. In the tables, the memory requirement is obtained by using the “size” command, a UNIX com- 
mand on an Ultra Sun, which gives the memory required to execute the code; the number of independent design vari- 
ables is equal to the number of nodes in an element. As the x-coordinates of the nodes are considered as the design 
variable, only a minimal increment in computer memory and time is observed in Table 1, for heat transfer elements. 
This is because, as a one-dimensional problem, the number of design variables is limited and this element does not 
involve coordinate transformation. 

In summary, the procedure presented here can effectively use ADIFOR to produce a code to calculate the deriva- 
tive terms required for the first-and second- order thermal sensitivity analyses. 


Table 1. Comparison of ADIFOR with finite difference for 1-D,/?- version element 


Element 

type 

Memory 

before 

ADIFORing 

(Bytes) 

Memory 

after 

ADIFORing 
for first 
derivatives 
(Bytes) 

Memory 

after 

ADIFORing 
for second 
derivatives 
(Bytes) 

CPU time 
for stiffness 
matrix 
evaluation 
of one 
element 
(sec) 

CPU time 
taken by 
ADIFOR for 
first 

derivatives 

(sec) 

CPU time 
taken by 
ADIFOR for 
second 
derivatives 
(sec) 

Number of 
independent 
variables 

1-D 

/?-version 

element 

13,706 

15,496 

20,864 

0.02 

0.03 

0.04 
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4, Numerical Studies 


The heat transfer problem studied here is one-dimensional and is modeled by p - version finite elements. It is used 
mainly to verify the accuracy of the sensitivity analyses, as it is simple and its exact derivatives can be analytically 
obtained. 


The steady-state heat transfer in one dimension is governed by the equation describing the conservation of the 
energy 



(ID 


where k is the thermal conductivity, T is the temperature, and q is the internal heat generation. In this study, the 
internal heat generation term is devised as 
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where the parameter a controls the steepness of the solution at x Q . As a result, the temperature distribution can be 
analytically solved by 


T(x) = (/-x)(tan l a(x - jc ) + tan ] (ax Q )) (13) 

which satisfies the boundary conditions, T( 0) = T(l) = 0. The value of /, a and x Q are specified as 1, 20 and 0.5, 
respectively. 

The analytical expressions of the first and second-order derivatives of T(x), with respect to the length of the 
domain, /, can then be derived as 
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where x 0 is replaced by 112, The change of domain, /, will induce changes in the locations of the mesh nodes and the 
location of the steepest temperature, x 0 . 

The p - version element is defined over a two-node, one-dimensional domain where the temperature distribution is 
interpolated as 


P 

t{ x) = 2> ( wr. 

i = 0 

Using the shape functions, N t of Demkowicz et al.: 10 

i = 0 
i = 1 

i = 2,4, 6,... 
i = 3,5,7, ... 
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where % is a mapping of the element coordinates onto a standard domain, E, e [-1, 11 . Note that the higher order 
shape functions have zero values at the end nodes. Hence, the coefficients of the linear shape functions, T 0 and T , . 
commonly used in an /z-version heat transfer element, do represent the temperatures at the nodes, while the coeffi- 
cients of higher order shape functions do not have physical meaning. It should also be noted that T(x) is a C° function 
throughout the temperature domain. 

The Galerkin method can be used here to obtain the matrix equation that approximates the solution of Eq. (11). 

KT = q (18) 

As described in Section 3, once a computer code is developed to solve the above equation, the ADIFOR-produced 
code can be used to solve Eqs. (9) and (10) for first- and the second-order derivatives, dT/db and d T/db , of the 
discrete temperature vector, T . 

To study accuracy, the heat transfer problem is discretized into 10 elements with 1 1 nodes. The total length, /, of 
the domain, is considered to be the design variable. The analytical shape derivatives of the temperature distribution 
with respect to the length of the domain can be found by Eqs.( 14)-(15), which are plotted in Fig. 2. On the other hand, 

the continuous expressions of the discrete temperature derivatives, dt/dl and d 2 t/dl 2 , can be interpolated by 


j , T '« ■ 1 11 M# 

i = o 

(19) 

and 


d 2 T £ d 1 T i 

dl dl 

(20) 

where is defined in Eqs. (7). The result shows in Fig. 2 that with the order of polynomials in N . 

greater than 3, the 

continuous expressions of the numerical derivatives are matched very well with those of analytical derivatives. 

The equations for the errors between the exact derivatives, dT/d/ and d 2 T/d/ 2 , and the /^-version finite element 

d 0 * 

solutions, -jrf and — -t , are given by 
dl dl 2 



(21) 

e(x) = £f(x) 

dl dl 

(22) 


To study the convergence of the finite element solution and derivatives with different meshes and approximation, the 
H l norms of the errors are introduced. For example, the H l norm of the error of the first order derivative is defined as 



N k 

where ||e||, k is an error in element k given by 
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similar definitions can be made for \\e\\ { and \\e\\ l . Three figurers are presented for errors of each type of solutions. 
The first two figures indicate the total H 1 errors of the solution when one of the parameters; the order of polynomials 
in approximation or the number of meshes, increases and the other remains unchanged. The third one is a picture of 
elemental errors in individual elements while the mesh size is being changed. 

Figures 4 to 6 show the convergence of the temperature solution. As discussed in Ref. 9, it is observed that con- 
vergence rate of p-version is faster than that of h-version finite element. Furthermore, Fig. 6, shows that elemental 
errors are not distributed uniformly. Particularly, high error is found near the neighborhood of high temperature gra- 
dient. 

Figures 7 to 9 show the convergences of the first-order temperature derivatives. The trend of the convergent rate 
of the temperature derivatives is similar to that of the temperature. Again, high error in temperature derivative is 
found near the neighborhood of high temperature gradient. However, it seems that the convergent rate of the temper- 
ature derivative is one order slower than that of the temperature solution and reducing mesh size may not improve the 
accuracy of the temperature derivative in some part of the domain, as shown in Fig. 9. 

Figures 10 to 12 show that neither reduction of mesh sizes nor increase of polynomial order can improve the con- 
vergence of the second-order temperature derivatives, measured by the H l error norm. This results are unexpected 
and require further study to fully understand such decay of convergence in high-order derivatives. 


5, Concluding Remarks 

The paper first proposes an implementation procedure to develop an ADIFOR-enhanced code for first- and sec- 
ond-order sensitivity analysis and then investigates the performance of the resultant code in terms of memory require- 
ment, computational accuracy and efficiency. The results have demonstrated that an user can follow the proposed 
procedure to develop an ADIFOR-enhanced code for calculating first-order and second-order derivatives without 
much knowledge of the original code and with only a minimal coding effort. The resultant code can calculate the first- 
order derivatives very efficiently. However, it needs a moderate increase of computer resources to support the calcula- 
tion of the second-order derivatives. The slow-down of the second-order sensitivity analysis is mainly caused by the 
computation of many second-order derivatives of elemental stiffness matrices, appearing on the right-hand side of the 
second-order sensitivity equation. The error analysis concludes that the convergence of the temperature derivatives 
show the same trend as that of the temperature solution. However, the convergence rate of the former is slower than 
that of the latter. Furthermore, in the cases studied, the h-adaptiveness fails to converge the temperature derivatives. 
Further research is needed to fully understand the error behavior of the temperature derivatives of which are calcu- 
lated by the p- version finite elements. 
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Fig. 4: Total H 1 error keeping p fixed and uniformly increasing h 


Fig. 5: Total H 1 error keeping h fixed and uniformly increasing p 
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Fig. 7: Total H l error keeping p fixed and uniformly increasing h 


Fig. 8: Total H l error keeping h fixed and uniformly increasing p 
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